// @HEADER
//***********************************************************************//
//    SiYuan: A numerical PDE solver                                     //
//    This Software is released under the BSD 2-Clause license detailed  //
//    in the file "LICENSE" in the top-level SiYuan directory            //
//***********************************************************************//
// @HEADER

#ifndef SIYUAN_CAHNHILLCHEMTERM_IMPL_HPP
#define SIYUAN_CAHNHILLCHEMTERM_IMPL_HPP

#include "Teuchos_TestForException.hpp"
#include "Phalanx_DataLayout.hpp"
#include "Sacado_ParameterRegistration.hpp"
//#include "Panzer_ExprEval_impl.hpp"

template<typename EvalT, typename Traits>
SiYuan::CahnHillChemTerm<EvalT, Traits>::
CahnHillChemTerm(Teuchos::ParameterList& p, const panzer::IntegrationRule&  ir)
{
  m_beta = p.get<double>("Beta");
  
  numQPs  = ir.dl_vector->extent_int(1);
  ir_degree = ir.order();

  m_rho = PHX::MDField<ScalarT, panzer::Cell, panzer::IP>("Concentration", ir.dl_scalar);
  this->addDependentField(m_rho);
  m_mu = PHX::MDField<ScalarT, panzer::Cell, panzer::IP>("Diffusion Potential", ir.dl_scalar);
  this->addDependentField(m_mu);

  m_freeEnergy = PHX::MDField<ScalarT, panzer::Cell, panzer::IP>("FreeEnergy", ir.dl_scalar);
  this->addEvaluatedField(m_freeEnergy);
  this->setName( "FreeEnergy" + PHX::print<EvalT>() );
}

// **********************************************************************
template<typename EvalT, typename Traits>
void SiYuan::CahnHillChemTerm<EvalT, Traits>::
postRegistrationSetup(typename Traits::SetupData sd,
                      PHX::FieldManager<Traits>& fm)
{
  this->utils.setFieldData(m_freeEnergy,fm);
  ir_index = panzer::getIntegrationRuleIndex(ir_degree,(*sd.worksets_)[0], this->wda);
}

// **********************************************************************
template<typename EvalT, typename Traits>
void SiYuan::CahnHillChemTerm<EvalT, Traits>::
evaluateFields(typename Traits::EvalData workset)
{
    for (std::size_t cell = 0; cell < workset.num_cells; ++cell) {
      for (std::size_t qp = 0; qp < numQPs; ++qp) {
      //  m_freeEnergy(cell, qp) = m_rho(cell, qp) * (m_rho(cell, qp)* m_rho(cell, qp)- m_beta*m_beta) 
      //    - m_mu(cell, qp);
            // FEnics:  100*c**2*(1-c)**2
        m_freeEnergy(cell, qp) = 200.0 * m_rho(cell, qp)* (1- m_rho(cell, qp)) * (1-2.0*m_rho(cell,qp))
          - m_mu(cell, qp);
      }
    }

}


// *****************************************************************************

#endif